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Abstract. We use the Bethe-Peierls method combined with the belief propagation 
algorithm to study the arctic curves in the six vertex model on a square lattice 
with domain-wall boundary conditions, and the six vertex model on a rectangular 
lattice with partial domain-wall boundary conditions. We show that this rather simple 
approximation yields results that are remarkably close to the exact ones when these are 
known, and allows one to estimate the location of the phase boundaries with relative 
little effort in cases in which exact results are not available. 
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1. Introduction 


Usually, boundary conditions have no effect on the equilibrium properties of statistical 
models in the thermodynamic limit. Still, boundary conditions may affect the 
equilibrium behaviour of strongly constrained models and, in particular, they may 
impose the macroscopic separation of phases. In such cases, the boundary between an 
internal spatial region, in which the configuration is typical of equilibrium with the bulk 
parameters, and an external spatial region, frozen by the boundary conditions is named 
the arctic curve. The external region is called arctic (for frozen) and the internal region is 
called temperate (for fluctuating). The arctic curve first appeared in the study of domino 
tilings of Aztec diamonds [H E] [S] [1] , then in lozenge tilings of large hexagons 016], 
and later in more general dimer [7] and vertex models 0 a El El dni EH [i2|. In these 
systems phase separation exists for a wide choice of fixed boundary conditions. 

The six vertex model was first introduced to model ferroelectricity [El Ell EH EEl 
EZj, and it has links with the models mentioned above. It is commonly defined on a 
square lattice with N x N vertices. (We will consider here the extension to rectangular 
systems with N x M vertices as well.) Arrows are placed along the links, with two 
possible orientations for each. The six vertex rule imposes that two arrows should point 
in and two should point out of each vertex. Energies and, consequently, statistical 
weights are assigned to each vertex. If complete arrow reversal symmetry is assumed, 
three parameters, a, b, c, characterise these weights, with the two first ones associated 
to local ferroelectric order and the last one linked to antiferroelectric order (see Fig. [^. 

The solution to the most general six vertex model with periodic boundary conditions, 
in the form of its bulk free-energy density, was given by Sutherland [E| who extended 
Lieb’s treatment with the Bethe Ansate technique to generic values of the vertex weights. 
The value taken by the following parameter (see Fig. 




+ 6 ^ — 
2ab 


( 1 . 1 ) 


distinguishes the three phases of the model: for A > 1 the equilibrium state has 
ferroelectric (FE) order, for A < —1 it has antiferroelectric (AE) order, and for |A| < 1 
it is disordered (D). The D-EE phase transition is first-order while the D-AF transition 
is characterized by an essential singularity reminiscent of the Kosterlitz-Thouless type. 

The six vertex model with domain wall boundary conditions was introduced in [E| 
to study correlation functions in exactly solvable models [20| . In a few words, domain- 
wall boundary conditions correspond to all arrows on the bottom and top boundaries 
being incoming while all arrows on the left and right boundaries being outgoing, see 
Fig.H (or vice versa). More recently, the interest in this model has been pushed by 
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its connection with algebraic combinatorics and the enumeration of alternating sign 
matrices. The partition function satisfies a recurrence relation [211 US], that leads to 
a determinant formula |2T], exploited by Korepin and P. Zinn-Justin (see also [22]) to 
derive exact expressions for the free-energy densities in the disordered and ferroelectric 
phases [23] • The antiferroelectric free-energy density was obtained in [23] by using a 
matrix model representation of the determinant. Although the phase diagram remains 
unaltered, and still given by the case with |A| = 1, the free-energy densities in the 
disordered and antiferromagnetic phases, are different from the ones found with periodic 
boundary conditions, even in the thermodynamic limit. Curiously enough, the free- 
energies take much simpler expressions as functions of the parameters a, b, c under the 
domain-wall boundary conditions. The D-FE transition becomes continuous while the 
D-AF remains of infinite order. 

The six vertex model can be mapped onto the problem of domino tilings on Aztec 
diamonds [IZllIEI for a particular choice of vertex weights m and, in consequence, it also 
admits an arctic curve separating an external finite density boundary with ferroelectric 
order from an internal region that is either disordered or antiferroelectrically ordered 
for |A| < 1 or A < — 1, respectively. On the free-fermion point (A = 0) with a = b the 
arctic curve is a circle Ell. On the less restrictive free-fermion cases with a ^ b the 
arctic curves are ellipses [E] . Numerical evidence for the existence of an arctic curve for 
general values of the parameters was presented in [251, [26|. Various analytic, though not 
yet proven to be exact (except for A = 0 where the derivation is rigorous), methods have 
been employed to determine the arctic curves in the disordered phase with A ^ 0 pun] 
and even in the antiferroelectric phase [ID¬ 
AS mentioned above, very powerful and interesting analytic methods allow one to 
understand the statics of the six vertex model, and the phase separation induced by 
special boundary conditions. A summary of modern methods applied to this problem 
can be found in m and a review on arctic curves in the six vertex model is given 
in na. Still, as soon as one lifts the integrability conditions satisfied by the six (and 
eight) vertex model, these techniques are no longer useful. Approximate methods, such 
as the Bethe-Peierls approximation |2E| and its modern versions, like the cavity method 
and the belief propagation algorithm [221 EDI ED E2] , can then be of great help to obtain 
the phaser diagram and equilibrium properties of generic vertex models. This method 
was applied in [331 El] to analyse the sixteen vertex model with parameters close to 
the ones of artificial spin-ice samples (331 EE] • Surprisingly enough, the method gave 
very accurate, sometimes even exact, results when applied to the integrable six and 
eight vertex cases. Moreover, cluster generalizations have been used to obtain phase 
diagrams and thermodynamical properties of vertex models with a much larger state 
space [3DEHlEnil3nil3DI32|, both in two and three dimensions. 

The work in [331 E3| concerned homogeneous systems (except for distinguishing 
between two different sublattices, when necessary) with periodic boundary conditions. 
Here we aim at extending this work to inhomogeneous systems, making it possible 
to describe the phase separation phenomenon induced by domain-wall boundary 
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conditions. This will imply an increase in the computational complexity by a factor 
N X M which, even in the large N and M limit, can be dealt with thanks to the 
belief propagation (BP) algorithm, which is particularly efficient at finding minima of 
a Bethe-Peierls free-energy even in inhomogeneous systems. As an example of results 
that we obtain in this work, for which there are no rigourous predictions, there is the 
behavior of the interface separating the disordered region from the AF phase in the case 
of parameters corresponding to the AF bulk phase. 

We end this introductory section by stating that domain wall boundary conditions 
could be easily imposed experimentally in artificial spin-ice samples I3S1ISE]- A possible 
way to do this would be to fabricate an array where the edge islands are different in 
some way (larger, or from a different magnetic material) in order to make them more 
stable. Modern visualization experimental methods |13l HU SSI SS] should then allow 
one to see the phase separating curves in the lab. 

The paper is organised as follows. In Sec. we recall the definition of the six vertex 
model and the properties that are of interest to our work. Section is devoted to a 
short description of the belief propagation method as applied to the six vertex model 
with domain wall boundary conditions. In Sec. we present our results. In Sec. we 
discuss several lines for future work. 

2. The six vertex model and its arctic curve 

2.1. Definition, phases and free-energies 

The six vertices defining the model together with their statistical weights are shown in 

Fig.[T} 


oji = a 0 J 2 = a UJ 3 = b UJ 4 = b uj^ = c ujq = c 

Figure 1: The six vertices in the six vertex model and their Boltzmann weights, oj\,... ,ojq. 

An alternative parametrisation of the Boltzmann weights is 

o = sin(A-|-r^) 6 = sin(A — r^) c = sin2ry (2.1) 

with A the ‘rapidity’ variable and rj the ‘crossing’ parameter. The parameter A = cos 2 ri 
serves to locate the phase transitions. 

In the disordered (D) phase |A| = | cos 2 / 7 1 < 1 and the parameters A and rj are 
constrained to p<X<tt — p and t) < p < n/2. In the free-fermion case A = 0 and 
r] = 7 r/ 4 . In the free-fermion and symmetric a = b case, A = 0, 77 = 7 r /4 and A = 7 r/ 2 . 
In the spin-ice case a = b = c implies A = 1/2, rj = tt/G and A = 7r/2. 
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Let us first focus on the square lattice model. The domain-wall boundary conditions 
(DWBC) are sketched in Fig.|^ The free-energy density per site in the disordered phase 
with these boundary conditions is [HI [23] 


fD _ 
J DW ~~ 


In 


a sin(A -|- rj) sin(A — rf) 
sin[Q;(A — rj)] 


with 


a = 


TT 

TT — 2ri 


( 2 . 2 ) 


When compared to the (much more involved) expression for periodic boundary 
conditions [TH] one finds > fp. The extension of this analysis to the 

antiferroelectric (AF) and ferroelectric (FE) phases shows that > fp^ in the AF 
phase Plj while fp^ = fp^ in the FE one [HI ES] • The macroscopic phase separation 
is intimately linked to the difference between the free-energies for different boundary 
conditions. The D-FE phase transition is continuous and the D-AF transition is of 
infinite order under DWBC. 



Figure 2: An example of domain wall boundary conditions. 


2.2. The arctic curve 

The arctic curve(s) delimiting the spatial regions with different ordering properties is 
(are) defined in the scaling limit in which the number of lines in each direction, on a 
NxN square lattice, tends to infinity while the lattice spacing, h, vanishes keeping the 
two linear size of the lattice fixed. 

In the case in which the bulk parameters select a disordered state, the phase 
separation induces five spatial regions in the sample: four external ferroelectric ones, 
close to the four corners and an internal one, that is the disordered region D. 

In the case in which the bulk parameters select the AF state, the phase separation 
occurs in two steps, in the sense that an intermediate disordered spatial region separates 
the external FE region fixed by the boundaries, and the internal AE region fixed by the 
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bulk parameters. There is therefore one arctic curve and one temperate (between D and 
AF regions) curve in this problem. 

Some analytic expressions for the interface between frozen and temperate regions 
in the six vertex model on a sqnare lattice are known and we snmmarise them below. 

2.2.1. Boundary between the FE and the D phase. Colomo and Pronko used boundary 
correlation functions to get a closed expression for the coordinates of the contact points 
between the arctic curve and the system’s boundaries [10] 

^ QCOt[Q(A - 7 ])] - COt(A + 7 ]) 

^ cot(A — rj) — cot(A + if) 

for generic |A| < 1. In particnlar, for A = 7 r /2 one finds k = 1/2 for all r^. 

The form of the arctic curve linking these points is highly non-trivial and for generic 
Tj non-algebraic. Some special cases have been known for more than a decade. 

Jockusch et al [1] proved the arctic circle theorem for the free-fermion case A = 0 
with a = 6 in which the arctic curve is just a circle, 



In the free fermion case, A = 0 and rj = vr/d, the arctic curve can be computed 
exactly for a 7 ^ 6 and it is an ellipse tangent to the four sides of the square 

(1 -S t^){x + y - 1)'^ + t^(l + = t^ (2.5) 

with 

t — - = tanh ~ • ( 2 - 6 ) 

For t = 1 one recovers the article circle. 

With the emptiness correlation function, and a conjecture on its behavior, Colomo 
and Pronko obtained an implicit expression for the arctic curve in the general A 7 ^ 
0 case P [IQI HZ]. In section 6.3 in [TO], eqs. (6.17)-(6.19) present a parametric 
representation of this curve. In general, the curve is a transcendental function. For 
special valnes of A, that correspond to what are called the roots-of-nnity in the quantum 
group context, the curve becomes algebraic nnnEj. In particular, at the spin-ice point 
[a = b = c and A = 1 / 2 ) the arctic curve is given by the portion of the ellipse 

{x -\- y — 2^ + 3(x — 1 /)^ — 3 = 0 (2.7) 

in which x,y E [0, 1 / 2 ]. The contact points in the lower left qnadrant are (1/2,0) and 
( 0 , 1 / 2 ), in other words, k = 1 / 2 . 

2.2.2. Boundary between the FE and AF phases. Colomo, Pronko and Zinn-Jnstin m 
used the emptiness correlation function method to estimate the arctic curve between the 
external frozen region and an intermediate disordered region that separates it from the 
bulk with AF order. As far as we know, there is no prediction for the internal boundary 
between the D and AF regions. 
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3. Method 


In this section we introduce the key ideas and equations of the belief propagation (BP) 
algorithm, which are needed to apply the method to the present problem. 

We label the lattice edges i,j, k,... and the vertices a,b,c,... (not to be confused 
with the vertex weight parameters used above). If an edge i is incident on a vertex 
a we say that the edge belongs to the vertex and we write i £ a. In the language of 
probabilistic graphical models (see e.g. [31] and refs, therein for the relationship between 
these models and statistical mechanics), edges (respectively, vertices) are variable nodes 
(resp. factor nodes) in a factor graph. In this notation the full Boltzmann weight can 
be written as 

i’ = Yli’aiSa), (3.1) 

a 

where G o} and s* represents the orientation of edge i, say, positive towards 

the right (up) for horizontal (vertical) links, as in |2S|- 

The BP algorithm is a message-passing algorithm for hnding the minima of an 
approximate (mean-field-like) variational free-energy (a generalization of the Bethe- 
Peierls one) which can be obtained by truncating the cumulant expansion of the entropy 
of the model, by keeping only edge and vertex contributions or, in other words, the free- 
energy of the cluster variational method (see [SI| for a recent review) with vertices as 
maximal clusters: 


O' — Ha(Sa)pa(Sa) + 

a Sa 


+ ksT 


EE Pa{Sa)\riPa{Sa) “ EE Pi{si) Inpi{si)\ , 

_ CL S d 1 S'i _ 


(3.2) 


where Ha{sa) = —kETlwipa^Sa) is the contribution of vertex a to the model Hamiltonian 
and Pa{sa) (respectively Pi{si)) is the probability distribution of vertex a (resp. edge i). 
The above variational free-energy has to be minimized with respect to the vertex and 
edge probability distributions, subject to the constraints of normalization 


J^Pa(Sa) = 1, Vo J^Pj(Si) = 1, VZ (3.3) 

Sa Si 


and marginalization 

'^PaiSa) = Pi{Si), Vo, VZ G O, (3.4) 

^a\i 


where Sa\i = {sj,j £ a, j i}. 

Minima of the variational free-energy correspond [HOI ED to fixed points of the BP 
algorithm. In the latter picture, a vertex a sends a message ma^i{si) to an edge z G o, 
which depends on the edge orientation. Up to a normalization constant, the edge and 
vertex probability distributions can be written as functions of the messages as 

b^a 

Pi{Si) OC JJma^i(Si), Pa{Sa) OC V'a('Sa) nn mb^i{si). (3.5) 

a3i i^a b3i 
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Rewriting the marginalization condition Eq. (3.4) in terms of messages one obtains, 
again up to a normalization constant, 

j¥^i 

ma^i{Si) OC ^'0a(Sa) JJ JJ (3-6) 

^a\i b3j 

the main iterative equation of the BP algorithm, the fixed point of which gives the 
message values at equilibrium, to be used in the computation of the edge and vertex 
probability distributions, Eq. (3.5), and of the free-energy, Eq. (3.2). In the following 
we shall often refer to the free-energy density / = J- /{NM). 

Domain-wall boundary conditions are easily represented in this scheme, by 
introducing boundary messages, sent by auxiliary vertices to boundary edges which 
vanish if and only if the edge orientation differs from the one imposed by the boundary 
conditions. In more detail, what we mean by this is the following. Consider a boundary 
edge i and assume that it is horizontal, on the right boundary of the lattice. This edge 
belongs to a single vertex a, which sits on its left side. Introduce the auxiliary vertex 
b, sitting on the right of the edge i. Since the polarization of i is constrained to be 
rightward (-|-1 in our representation), we set 1) = 0. 


4. Results 

In this Section we present the results that we obtained by using the BP algorithm 
explained in Sec. 

4 -1. Square lattice 

We first study the model on a square lattice for which analytic expressions for the 
interfaces are known. 


4.TT The free-energy. In |3S] the Bethe-Peierls approximation was used to study 
the equilibrium properties of the sixteen and, in particular, the six vertex model with 
periodic boundary conditions. Defining the six vertex model on a tree, the exact EE free- 
energy density was found, while in the D and AF phases very good approximations to 
the exact expressions, that get better and better far from the transition, were derived. 
Although the free-energy densities obtained with the Bethe-Peierls approximation in 
the infinite size limit (that in the homogeneous case can be treated analytically) are not 
exact, the location of the different transition lines (between D and EE, and D and AF, 
phases) are. 

First, we check that the BP free-energy densities of the six-vertex model with 
DWBC satisfies fow > fp for a square lattice system. In Fig. we show the approach 
to the thermodynamic limit of the free-energy density at the spin-ice point o = 6 = c = 1 
(A = 1/2). The numerical data are fitted with the form 

fSwi^) — fSwi^ ^ • 


(4.1) 
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We find ~ 1.6 — 1.7 and fEwi^ — —0.246. The estimate for the free- 

energy density in the thermodynamic limit is larger than the one for periodic boundary 
conditions, fp = — ln(3/2) ~ —0.404, that is Pauling’s result also found on the tree. 
We note that the free-energy-density for domain-wall boundary conditions on the tree 
is slightly larger than the exact result, fEw = — ln(3-s/3/4) ~ —0.262 j2S]- The fact 
that the BP approximation slightly over-estimates the free-energy (with respect to the 
exact result) was also found for periodic boundary conditions [HHj . 



Figure 3: BP free-energy density vs. lattice size N for the six vertex model on a square lattice 
with domain-wall boundary conditions at the spin-ice point, a = b = c = 1 (data points joined 
by straight lines). The asymptotic estimate, fEw — —0-246, is shown with a dotted horizontal 
lines. The exact value, fEw — “0-262, lies beyond the vertical range of variation of this graph. 

In Fig. 1^ we compare our free-energy density with the exact one [23], in the D and 
FE phases. In particular we consider a square lattice with N = 32, we set b = c = 1 and 
vary o in the range [1,3]. Notice that a — 1 corresponds to the spin-ice point A = 1/2, 
while at a = 2 (A = 1) the system is expected to undergo (in the thermodynamical 
limit) a continuous phase transition between the D and FE phases. Our finite size 
analysis shows instead a continuous variation from a D bulk behaviour for a < 2 to an 
EE bulk behaviour for a > 2. In agreement with data reported in Eig. the corrections 
to our free-energy density that we obtain upon increasing N from 32 to 1000 are of 
order 10“^, hence not resolvable on the scale of Eig. [ij 

The exact FE free-energy of the two-dimensional model in the A —>■ cxo limit is 
fEw — fp^ = — Ino and it is shown with a blue line in Fig. The BP algorithm 
should also give this result in the infinite size limit. We ascribe the deviation from of 
the black curve in Fig. |^from the asymptotic limit to finite size effects. 

A similar analysis is shown in Fig. in the D and AF phases (for the exact free- 
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Figure 4: Flee-energy density vs. a for the six vertex model on a square lattice with parameters 
b = c = 1 and under domain-wall boundary conditions. The BP data for a system with linear 
size N = 32 (black line) should be compared to the exact result [231 thermodynamical 

limit (red line: D phase, blue line: FE phase). 


energy density in the AF phase see [21] )• We set N = 32, a — b = 1 and vary c in 
the range [1,3]. Again c = 1 corresponds to the spin-ice point A = 1/2, while at 
c = 2 (A = — 1) the system is expected to undergo (in the thermodynamical limit) a 
continuous phase transition between the D and AF phases. Here our finite size results 
show a reminiscence of the infinite size phase transition. For c G [2.08, 2.13] we find both 
a (stable) AF phase and a (metastable) D phase (the difference in free-energy density 
between these two phases is ~ 10“^, not resolvable on the scale of Fig. [^, while for 
c < 2.08 we have only the D phase and for c > 2.13 only the AF phase is found. This 
allows us to locate a finite size phase transition at c ~ 2.08. Upon increasing N the 
transition point seems to tend towards c = 2, but it is hard to prove that the transition 
line is captured exactly as N ^ oo. 

4 . 1 . 2 . Arctic circle. We work out the case A = 0 and o = 6 in which the arctic curve 
should be a circle. In Fig. |^we plot the polarization (sj) of the horizontal edges for a 
1024 X 1024 lattice, with (following the convention introduced in Sec. above) s, = +1 
(respectively, —1) for a rightward (resp. leftward) edge. The exactly known arctic circle 
is drawn as a white solid line. The agreement is remarkable, considering that we are 
applying a mean-field like technique to a system in its critical state. The temperate 
region is slightly overestimated by the BP approximation, though this might be just 
a finite size effect. This is better observed in Fig. [^ where we plot the polarization 
(from now on, the polarization of the horizontal edge at position {x,y) will be denoted 
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Figure 5: Free-energy density vs. c for the six vertex model on a square lattice with parameters 
a = b = 1 and under domain-wall boundary conditions. The BP data for a system with linear 
size N = 32 (black line) should be compared to the exact result |23l 124) in the thermodynamical 
limit (red line: D phase, blue line: AF phase). 


by Px{x,y)) of the horizontal edges at a fixed vertical coordinate y = 1/8, for different 
lattice sizes, together with the left intersection of the exact arctic circle with y = 1/8. 

Given that the system is critical, we expect simple power-laws for the scaling of 
its observables. We denote by XQ{y) the (exactly known) arctic curve. Observing that 
in the frozen regions on the right of Fig. the polarization Px{x,y) tends to 1, we 
have considered the deviation 1 — px{x,y) as a function of A^“(x — XQ{y)) for fixed y 
and a: > 1/2 and looked for data collapse. The best results have been obtained with 
a ~ 1/6 and the corresponding plot is reported in Fig. with the vertical axis both in 
linear (left panel) and log (right panel) scale. A good collapse is obtained for x > Xo{y) 
suggesting that our estimate of the arctic curve is very close to the exact one and that 
the scaling law holds in the frozen phase only. 

4 . 1 . 3 . Arctic ellipses. We now work out the case A = 1/2 and o = 6 = c in which the 
arctic curve should be an ellipse. In Fig. |^we plot, as in the previous subsection, the 
polarization of the horizontal edges and the exactly known arctic ellipse. The agreement 
is still remarkable. Our numerical data are consistent with the exact k = 1/2 values 
for the contact points. However, also in this case the temperate region seems slightly 
overestimated by the BP approximation. Indeed, in Fig. [T^ we display the polarisation 
of the horizontal edges as a function of the horizontal coordinate for a cut at constant 
vertical coordinate, y = 1/8. The right panel shows a zoom close to the x value where 
the exact article curve lies. One sees that the finite N curves obtained with the BP 
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0 0.2 0.4 0.6 0.8 1 


Figure 6: BP polarization of the horizontal edges for the free-fermion case, A = 0, with a = 6, 
on a square lattice with 1024 x 1024 lattice vertices. The white line is the exact arctic circle. 



Figure 7: BP polarization of the horizontal edges for the free-fermion case, A = 0, with a — b^ 
for ^ = 1/8 and different system sizes N. The vertical line denotes the left intersection of the 
exact arctic circle with ^ = 1/8. 


algorithm suggest a location of the sharp boundary at a value of x that will be closer 
to the edge of the system than the exact one. 









P/x,y) 
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Figure 8: BP polarization (deviation from the value in the frozen phase) of the horizontal 
edges vs. N^{x — xo{y)) for the free-fermion case, A = 0, with a — b, for y = 3/4, a — 1/6 
and different system sizes N. Vertical axis in linear (left panel) and log (right panel) scale. 



Figure 9: BP polarization of the horizontal edges at the spin-ice point, a = 6 = c, A = 1/2. 


4.1.4- Double interfaces in the AF phase. For bulk parameters such that the system is 
in the AF phase, the domain wall boundary conditions impose a frozen external region, 
an intermediate disordered ‘buffer’ and an internal region with AF order. This is shown 
in Fig. 11 for different values of c > 2 and o = 6 = 1 (the phase transition between the 


bulk D and AF phases is located at c = 2). 

The shrinking of the internal AF region as c decreases towards 2 suggests a 




















P^(x,y) 
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Figure 10: Left panel: BP polarization of the horizontal edges at the spin-ice point, a = b = 

A = 1/2, for y = 1/8 and different system sizes N. The vertical line denotes the left 
intersection of the exact arctic ellipse with y = 1/8. Right panel: magnification close to 
the arctic ellipse. 


continuous transition: indeed, this transition is known to be continuous, and in 
particular of infinite order [24] (a feature which of course cannot be reproduced by 
our mean-field-like approach). 



0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 


Figure 11: BP polarization of the horizontal edges for A < — 1 with a = b = 1 and 
c = 3, 2.5, 2.05 from left to right. The system size is A = 256. 

While the variation of the AF-D boundary with c is evident from Fig. [m it is not 
the same for the D-frozen boundary (the arctic curve). We can however understand its 
trend from Fig. [T^ where we have plotted the polarization of the horizontal edges at a 
fixed height y = 1/8 for various values of c. Looking at the points where the polarization 
saturates it is now evident that the arctic curve gets larger as c decreases, both for c > 2 
(bulk AF phase, curves characterized by oscillations in the center) and for c < 2 (bulk 
D phase, smooth curves). 
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Figure 12: BP polarization of the horizontal edges at a = 6 = 1 meaning A = (2 — c^)/2, 
system size N — 256, and y — 1/8, for different c values. For c = 3, 2.5, 2.05, A < —1 and the 
system prefers the AF configuration in the bulk. For c=l, 0.1,0<A<1 and the system 
prefers the D configuration in the bulk. 


4-2. Rectangular systems with partial domain wall boundary conditions 

In the framework of the BP algorithm it is straightforward to consider rectangular 
lattices and different types of boundary conditions. In order to illustrate this, in the 
present Section we discuss the six vertex model on rectangular lattices, with partial 
domain wall boundary conditions (pDWBC). These are a generalization of DWBC to 
rectangular systems [HI 09]. Consider a rectangular lattice with N x M vertices (in 
the following we shall assume A > M to fix ideas) and impose DWBC-like boundary 
conditions on three sides of the lattice, leaving the 4**^ side free. Following the convention 
we have set up in Fig. we will have inward-pointing edges on the bottom side and 
outward-pointing edges on the left and right sides. The ice rule then implies that on the 
top side we must have M inward-pointing edges and N — M outward-pointing edges 
(hence for M = N DWBC are recovered). The condition that the top side is left free is 
equivalent to sum over all possible distributions of its M inward-pointing edges. 

The effect of pDWBC is illustrated in Fig. [^with a 48 x 32 system and different 
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values of the parameters, corresponding to the FE, D and AF phases from left to right. 
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Figure 13: BP polarization of the horizontal edges for a 48 x 32 rectangular system and 
{a,b,c) = (3,1,1) (FE phase, left panel), (1,1,1) (D phase, center panel) and (1,1,2.5) (AF 
phase, right panel). 

In the FE phase the effect of pDWBC is trivial: the interface between the two 
symmetric homogeneous phases remains at 45° with respect to the lattice boundary, and 
the M inward-pointing edges of the top side are, in the most probable arrangement, 
placed on the right. In the D and AF phase the phase separation phenomena and the 
arctic line are still observed, although with a different geometry which, qualitatively, 
can be thought of as obtained by cutting the top 16 rows of a square 48 x 48 lattice. 

The most interesting phenomena are observed by varying the aspect ratio in the 
AF phase, as shown in Fig. 

As the aspect ratio decreases the inner AF region is first split in 2 regions, which 
then get farther from each other and eventually disappear. 

5. Outlook 

We studied the spatial organisation of vertices in the six vertex model with domain wall 
boundary conditions by using an approximate method, the belief propagation or Bethe- 
Peierls technique. We found that, although the method is not expected to give exact 
results for the location of the various arctic and temperate curves, the forms obtained 
are remarkably close to the analytic ones, when these are known. The advantage of this 
method is that it is quite simple to implement and that it can be applied to situations 
in which the various exact methods used in the literature do not necessarily apply. 
For instance, any other type of fixed boundary conditions [HO], or asymmetric weights 
{ui 7 ^ UJ 2 , 0 i >3 7 ^ a; 4 , (Us 7 ^ uq) could be easily dealt with using the same technique. 

Our study has been static. It would be interesting to analyse the dynamics of 
formation of such structures allowing for the presence of a small density of defects, 
e.g. along the lines in ISD E2l ESI EU, and with protocols as the ones used 
experimentally [SSI ESI E7] . The analysis of finite size systems with various aspect 
ratios will be especially relevant to compare with experiments. 

Artificial spin-ice samples can be imaged in real space with various magnetic 
microscopy such as Lorentz microscopy |13] , photoelectron emission microscopy gmis] 
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Figure 14: BP polarization of the horizontal edges in the AF phase for a = 6 = 1, c = 2.5 
for rectangular systems of size 64 x 32 (left panel), 96 x 32 (center panel) and 128 x 32 (right 
panel). 


and magnetic force microscopy [46]. With these techniques one could access the spatial 
structures that we discussed in this work. 
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